Generalizing the wavelet-based multifractal formalism to vector-valued random fields: 
application to turbulent velocity and vorticity 3D numerical data 



Pierre Kestener^'^ and Alain Arneodo^ 

^CEA-Saclay, DSM/DAPNIA/SEDI, 91191 Gif-sur-Yvette, France 
^Laboratoire de Physique, Ecole Normale Superieure de Lyon, 46 allee d'ltalie, 69364 Lyon cedex 07, France 

(Dated: February 2, 2008) 

We use singular value decomposition techniques to generalize the wavelet transform modulus 
maxima method to the multifractal analysis of vector-valued random fields. The method is calibrated 
on synthetic multifractal 2D vector measures and monofractal 3D fractional Brownian vector fields. 
We report the results of some application to the velocity and vorticity fields issued from 3D isotropic 
turbulence simulations. This study reveals the existence of an intimate relationship between the 
singularity spectra of these two vector fields which are found significantly more intermittent than 
previously estimated from longitudinal and transverse velocity increment statistics. 
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The multifractal formalism was introduced in the con- 
text of fully-developed turbulence data analysis and mod- 
eling to account for the experimental observation of some 
deviation to Kolmogorov theory (K41) of homogenous 
and isotropic turbulence The predictions of vari- 
ous multiplicative cascade models, including the weighted 
curdling (binomial) model proposed by Mandelbrot 0, 
were tested using box-counting (BC) estimates of the 
so-called f{a) singularity spectrum of the dissipation 
field Alternatively, the intermittent nature of the 
velocity fluctuations were investigated via the computa- 
tion of the D{h) singularity spectrum using the structure 
function (SF) method 0. Unfortunately, both types of 
studies suffered from severe insufficiencies. On the one 
hand, they were mostly limited by one point probe mea- 
surements to the analysis of one (longitudinal) velocity 
component and to some ID surrogate approximation of 
the dissipation |5|. On the other hand, both the BC 
and SF methodologies have intrinsic limitations and fail 
to fully characterize the corresponding singularity spec- 
trum since only the strongest singularities are a priori 
amenable to these techniques 0. In the early nineties, 
a wavelet-based statistical approach was proposed as a 
unified multifractal description of singular measures and 
multi-afRne functions |^. Applications of the so-called 
wavelet transform modulus maxima (WTMM) method 
have already provided insight into a wide variety of prob- 
lems, e.g., fully developed turbulence, econophysics, me- 
teorology, physiology and DNA sequences [Tj, |3 . Later 
on, the WTMM method was generalized to 2D for mul- 
tifractal analysis of rough surfaces , with very promis- 
ing results in the context of the geophysical study of 
the intermittent nature of satellite images of the cloud 
structure 0, 0| and the medical assist in the diag- 
nosis in digitized mammograms [HI . Recently the 
WTMM method has been further extended to 3D analy- 
sis and applied to dissipation and enstrophy 3D numeri- 
cal data issue from isotrorac turbulence direct numerical 
simulations (DNS) [H [l3|. Thus far, the multifractal 



description has been mainly devoted to scalar measures 
and functions. In the spirit of a preliminary theoretical 
study of self-similar vector-valued measures by Falconer 
and O'Neil 0|, our objective here is to generalize the 
WTMM method to vector-valued random fields with the 
specific goal to achieve a comparative 3D vectorial mul- 
tifractal analysis of DNS velocity and vorticity fields. 

Let us note V(x = {xi,X2, ..,Xd)), a vector field with 
square integrable scalar components Vj{x)^j = 1,2, 
Along the Hue of the 3D WTMM method Hi Q, let us 
define d wavelets tpi{^) = d(j){x)/dxi for i = 1,2, .., d re- 
spectively, where i^(x) is a scalar smoothing function well 
localized around |x| — 0. The wavelet transform (WT) 
of V at point b and scale a is the following tensor [l^ : 



T^[V](b,a) 



, (1) 



where 



r^jy,](b,a) = a-'' / d''r^,{a~\r~h))Vj{v). (2) 



In order to characterize the local Holder regularity of V, 
one needs to find the direction that locally corresponds to 
the maximum amplitude variation of V. This can be ob- 
tained from the singular value decomposition (SVD) [l^ 
of the matrix {T^A^A) (Eq. O): 



(3) 



where G and H are orthogonal matrices {G'^G — H-^H = 
Id) and S = diag{<Ti, a^, .., Od) with Ui > 0, for I < i < d. 
The columns of G and H are referred to as the left 
and right singular vectors, and the singular values of 
[V] arc the non-negative square roots at of the d 
eigenvalues of T.0[V]^T.0[V]. Note that this decompo- 
sition is unique, up to some permutation of the cJi's. 



2 



The direction of the largest ampHtude variation of V, 
at point b and scale a, is thus given by the eigenvec- 
tor Gp(b, a) associated to the spectral radius p{h,a) = 
maxj (Jj{h, a). One is thus led to the analysis of the vec- 
tor field T^^p[V](b,a) = p(h, a)G n(h a) Following the 
WTMM analysis of scalar fields |l[i3[lj, let us define, 
at a given scale a, the WTMM as the position b where 
the modulus A^^[V](b, a) = |T^,p[V](b, a)| = p(b, a) is 
locally maximum along the direction of Gp(b, a). These 
WTMM lie on connected {d — 1) hypersurfaces called 
maxima hypersurfaces (see Figs l^b and[2t). In theory, 
at each scale a, one only needs to record the position of 
the local maxima of (WTMMM) along the maxima 
hypersurfaces together with the value of A^^[V] and the 
direction of Gp. These WTMMM are disposed along con- 
nected curves across scales called maxima lines living in 
a (d-l- 1) space (x, a). The WT skeleton is then defined as 
the set of maxima lines that converge to the (xi, X2, Xd) 
hyperplane in the limit a 0+ (see Fig. EJi). The local 
Holder regularity of V is estimated from the power-law 
behavior 7W^[V](£ro (a)) ~ along the maxima line 

£ro(fl) pointing to the point rp in the limit a — > 0"*", 
provided the Holder exponent /i(ro) be smaller than the 
number of zero moments of the analyzing wavelet 
HI . As for scalar fields HUES . the tensorial WTMM 
method consists in defining the partition functions: 



Z{q,a) 



E 

CeC{a) 



(A^^[V](r,a))* ^ a 



(4) 



where q G M and C{a) is the set of maxima lines that 
exist at scale a in the WT skeleton. Then by Legen- 
dre transforming t((;), one gets the singularity spectrum 
D{h) = miiiq{qh — T{q)), defined as the Hausdorff dimen- 
sion of the set of points r where h{r) = h. Alternatively, 
one can compute the mean quantities: 

h{q,a)= \n\M^[V]{r,a)\ W^[V]{q,£,a) , 

£e£(a) (5) 
D{q,a)^ J2 W^[V]{q,C,a) \n{W^[V]{q,£,a)) , 

CeC(a) 

where W^[Y]{q,C,a) = {M^[Y]{Y,a)Y / Z{q,a) is a 
Boltzmann weight computed from the WT skeleton. 
From the scaling behavior of these quantities, one 
can extract h(q) = \\mg^^Q+ h{q,a) / hia and D{q) = 
lim(j^Q+ D{q,a) / \na and therefore the D{h) spectrum. 

As a test application of this extension of the WTMM 
method to the vector situation, let us consider the self- 
similar 2D vector measures supported by the unit square 
defined in Ref. . As sketched in Fig. ^ from step n 
to step n -|- 1, each square is divided into 4 identical sub- 
squares and for each of these sub-squares, one defines a 
similitude Si that transforms the vector V^") at step n 
into the vector V^-"^^^ for the sub-square i at step n -f 1. 
The (T-additivity property of positive scalar measures is 



(a) 



(b) 



(c) 



FIG. 1: First construction steps of a singular vector- valued 
measure supported by the unit square. The norm of the four 
similitude Si axe pi = p4 — 1/2, p2 — 2 and ps = 1 ,15]. 



now replaced by the vectorial additivity condition V*^") = 
J2i=i^'i^^^^ ■ A straightforward calculation yields the 
following analytical expression for the partition function 
scaling exponents r(g) (Eq. Q): 



T{q) = -log2{pl+pl+pl+pl) 



(6) 



where pi {i ^ I to 4), are the norms of the similitudes 
Si. Note that this formula is identical to the theoret- 
ical spectrum of a non-conservative scalar multinomial 
measure distributed on the unit square with the weights 
Pi EE Q| ■ Indeed, if the construction process in Fig. 
is conservative from a vectorial point of view, it does not 
conserve the norm of the measure: = 4 > 1. 

From Legendre transforming Eq. ©, one gets a D{h) 
singularity spectrum with a characteristic multifractal 
single- humped shape (see Fig.|3i) supported by the inter- 
val [/imin, /imax] = [- 1 - loga (max^ ) , - 1 - loga (min^ )] 
and whose maximum Dp — — t(0) = 2 is the signature 
that the considered vector-valued measure is almost ev- 
erywhere singular on the unit square. 

In Fig. |2] are illustrated the main steps of our tensorial 
WT methodology when applied to 16 (1024)^ realizations 
of a random generalization of the vectorial multiplicative 
construction process described in Fig. Q Focusing on the 
central (128)^ sub-square, we show the singular vector- 
valued measure (Fig. EJl) and the corresponding WTMM 
chains computed with a first order analyzing wavelet at 
two different scales (Figs and|2t). On these max- 
ima chains, the black dots correspond to the location 
of the WTMMM at these scales. The size of the ar- 
rows that originate from each black dot is proportional 
to the spectral radius /o(b, a) and its direction is along 
the eigenvector Gp(b,a). When linking these WTMMM 
across scales, one gets the set of maxima lines shown 
in Fig. as defining the WT skeleton. In Fig. 13 are 
reported the results of the computation of the multifrac- 
tal spectra (annealed averaging). As shown in Fig. 
Z{q,a) (Eq. |@J) display nice scaHng behavior over four 
octaves (when plotted versus a in a logarithmic represen- 
tation), for 9 e] — 2,4[ for which statistical convergence 
turns out to be achieved. A linear regression fit of the 
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FIG. 2: 2D WT analysis of the 2D vector- valued self-similar 
measure shown in Fig. but with systematic random permu- 
tation of the Si. V is a first-order analyzing wavelet (0(r) is 
the Gaussian), (a) 32 grey-scale coding of the central (128)^ 
portion of the original (1024)^ field. In (b) o = 2^aw and 
(c) a — 2'^ aw, are shown the maxima chains; from the local 
maxima (WTMMM) of A4^ along these chains (■) originates 
a black arrow whose length is proportional to and direc- 
tion is along T^,p[V]. (d) WT skeleton obtained by linking 
the WTMMM across scales, aw ~ 7 (pixels) is the charac- 
teristic size of ip at the smallest resolved scale. 



data yields the nonlinear T{q) spectrum shown in Fig.|3|;, 
in remarkable agreement with the theoretical spectrum 
(Eq. 10). This multifractal diagnosis is confirmed in 
Fig. ^ where the slope of h{q, a) (Eq. Q) versus logj a, 
clearly depends on q. From the estimate of h(q) and 
D{q) (Eq. (jSJ), one gets the single-humped D{h) curve 
shown in Fig. which matches perfectly the theoretical 
D{h) spectrum. In Fig. O we have reported for com- 
parison, the results obtained when using a box-counting 
(BC) algorithm adapted to the multifractal analysis of 
singular vector- valued measures {14. ,15. _18J . There is no 
doubt that BC provides much poorer results, especially 
concerning the estimates of r(g), h{q) and D{q) for nega- 
tive q values. This deficiency mainly results from the fact 
that the vectorial resultant may be very small whereas 
the norms of the vector measures in the sub-boxes are 
not small at all. The results reported in Fig. |2| bring 
the demonstration that our tensorial WTMM methodol- 
ogy paves the way from multifractal analysis of singular 
scalar measures to singular vector measures. 

In Fig. 21 are reported the results of the application 
of our tensorial WTMM method to isotropic turbulence 
DNS data obtained by Leveque. This comparative 3D 
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FIG. 3: Multifractal analysis of the 2D vector- valued random 
measure field using the 2D tensorial WTMM method (□) and 
BC techniques (■). (a) logj .Z(g, a) vs logj a; (b) h{q,a) vs 
logj a; the solid lines correspond to linear regression fits over 
aw ^ a, < 2^aw ■ (c) r(g) vs q\ the solid line corresponds to 
the theoretical prediction (Eq. Q). (d) D(li) vs h\ the solid 
line is the Legendre transform of Eq. @. 

multifractal analysis of the velocity (v) and vorticity [ui) 
fields corresponds to some averaging over 18 snapshots 
of (256)3 DNS run at Rx = 140. As shown in Figs. Hi 
and^la, both the Z{q,a) and h(q,a) partition functions 
display rather nice scaling properties for q = —4 to 6, 
except at small scales (a < where some cur- 

vature is observed in the log-log plots likely induced 
by dissipation effects P, Linear regression fit of 

the data (Fig. 0^) in the range < a < 

yields the nonlinear Tv{q) and r^^iq) spectra shown in 
Fig. the hallmark of multifractality. For the vortic- 
ity field, Ti^{q) is a decreasing function similar to the 
one obtained for the singular vector-valued measure in 
Fig. Ot; hence h{q){— dT{q)/dq)< and the support 
of the D{h) singularity spectrum expands over negative 
h values as shown in Fig. In contrast Tv(g) is an 
increasing function which implies that h{q) > as the 
signature that v is a continuous function. Let us point 
out that the so-obtained Tv (q) curve significantly departs 
from the linear behavior obtained for 18 (256)'^ realiza- 
tions of vector-valued fractional Brownian motions B^/^ 
of index H = 1/3, in good agreement with the theo- 
retical spectrum T-Qi/3{q) = q/3 — 3. But even more 
remarkable, the results reported in Fig. for h{q, a) 
suggest, up to statistical uncertainty, the validity of the 
relationship h^{q) = h^{q) — 1- Actually, as shown in 
Fig. Eli, Di^{h) and D^{h) curves are likely to coincide 
after translating the later by one unit on the left. This 
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FIG. 4: Multifractal analysis of Leveque DNS velocity (•) 
and vorticity (o) fields (d — 3, 18 snapshots) using the ten- 
sorial 3D WTMM method; the symbols (■) correspond to a 
similar analysis of vector-valued fractional Brownian motions, 
■gH-i/3 ^^-j logj -E((jr, a) vs logj a; (b) hui{q,a) vs log 2 a and 
hv[q,a) — logj a vs logj a; the solid and dashed lines corre- 
spond to linear regression fits over 2^'^ aw ^ a < 2'^'^ aw- (c) 
Tv(q), T^{q) and r^i/^iq) vs g; (d) D^ih + 1), Du,{h) vs h; 
the dashed lines correspond to log-normal regression fits with 
the parameter values CJ ~ 0.049 and C" = 0.055; the dotted 
line is the experimental singularity spectrum (C2 ~ 0.025) 
for ID longitudinal velocity increments ||T8|. 



is to our knowledge the first numerical evidence that the 
singularity spectra of v and u) might be so intimately re- 
lated: Dv{h + 1) — Di^{h) (a result that could have been 
guessed intuitively by noticing that u; = V A v involves 
first order derivatives only). Finally, let us note that, for 
both fields, the T{q) and D{h) data are quite well fitted 
by log-normal parabolic spectra 19]: 



T{q) = -Co - Ciq - C29V2 , 
D(/i)= Co-(/i + Ci)V2C2. 



(7) 



Both fields are found singular almost everywhere: Cg = 
-rv(g = 0) = i:>v((7 = 0) = 3.02 ± 0.02 and = 
3.01 ± 0.02. The most frequent Holder exponent h^q = 
0) — —C\ (corresponding to the maximum of D{h)) 
takes the value -C^ ~ -Cf + 1 = 0.34 ± 0.02. In- 
deed, this estimate is much closer to the K41 prediction 
h = 1/3^ than previous experimental measurements 
{h — 0.39 ± 0.02) based on the analysis of longitudi- 
nal velocity fluctuations Consistent estimates are 
obtained for C2 (that characterizes the width of D{h)): 
CJ = 0.049 ± 0.003 and = 0.055 ± 0.004. Note that 
these values are much larger than the experimental es- 
timate C2 = 0.025 ± 0.003 derived for ID longitudinal 



velocity increment statistics |l9j . Actually they are com- 
parable to the value C2 = 0.040 extracted from experi- 
mental transverse velocity increments [19b]. 

To conclude, we have generalized the WTMM method 
to vector-valued random fields. Preliminary applications 
to DNS turbulence data have revealed the existence of 
an intimate relationship between the velocity and vortic- 
ity 3D statistics that turn out to be significantly more 
intermittent than previously estimated from ID longitu- 
dinal velocity increments statistics. This new method- 
ology looks very promising to many extents. Thanks to 
the SVD, one can focus on fluctuations that are locally 
confined in 2D (min^ ai = 0) or in ID (the two smallest 
(Ti are zero) and then simultaneously proceed to a mul- 
tifractal and structural analysis of turbulent fiows. The 
investigation along this line of vorticity sheets and vortic- 
ity filaments in DNS is in current progress. We are very 
grateful to E. Leveque for allowing us to have access to 
his DNS data and to the CNRS under GDR turbulence. 



[7] 
[8] 
[9] 
[10] 

[11] 

[12] 
[13] 

[14] 
[15] 

[16] 

[17] 



U. Frisch, Turbulence (Cambridge Univ. Press, Cam- 
bridge, 1995). 

B. B. Mandelbrot, J. Fluid Mech. 62, 331 (1974). 

C. Meneveau and K. R. Sreenivasan, J. Fluid Mech. 224, 
429 (1991). 

G. Parisi and U. Frisch, in Turbulence and Predictabil- 
ity in Geophysical Fluid Dynamics and Climate Dynam- 
ics, edited by M. Ghil et al. (North-Holland, Amsterdam, 
1985), p. 84. 

E. AureU et al, J. Fluid Mech. 238, 467 (1992). 
J. F. Muzy, E. Bacry, and A. Arneodo, Phys. Rev. E 
47, 875 (1993); Int. J. of Bifurcation and Chaos 4, 245 
(1994); A. Arneodo, E. Bacry, and J. F. Muzy, Physica 
A 213, 232 (1995). 

A. Arneodo et al., Ondelettes, Multifractales et Turbu- 
lences : de I'ADN aux croissances cristallines (Diderot 
Editeur, Art et Sciences, Paris, 1995). 
The Science of Disasters : climate disruptions, heart at- 
tacks and market crashes, edited by A. Bunde, J. Kropp, 
and H. Schellnhuber (Springer Verlag, Berlin, 2002). 
A. Arneodo, N. Decoster and S. G. Roux, Eur. Phys. 
J. B 15, 567 (2000); N. Decoster, S. G. Roux, and A. 
Arneodo, Eur. Phys. J. B 15, 739 (2000). 
A. Arneodo, N. Decoster, and S. G. Roux, Phys. Rev. 
Lett. 83, 1255 (1999). S. G. Roux, A. Arneodo, and N. 
Decoster, Eur. Phys. J. B 15, 765 (2000). 
A. Arneodo et al.. Advances in Imaging and Electron 
Physics, 126, 1 (2003). 

P. Kestener et al, Image Anal. Stereol. 20, 169 (2001). 
P. Kestener and A. Arneodo, Phys. Rev. Lett. 91, 194501 
(2003). 

P. Kestener, Ph.D. thesis, University of Bordeaux I, 2003. 
K. J. Falconer and T. C. O'Neil, Proc. R. Soc. Lond. A 
452, 1433 (1996). 

G. H. Golub and C. V. Loan, Matrix Computations, 2nd 
ed. (John Hopkins University Press, Baltimore, 1989). 
Note that if /li(ro) are the Holder exponents of the d 
scalar components Vi(ro) of V, then /i(ro) = min^ hiijco). 



5 



[18] We recall the relationship h = a — d, between the Holder 
exponent h and the singularity exponent a generally ob- 
tained by BC techniques {D(a — d) — /(a)). 

[19] (a) A. Arneodo, S. Manneville, and J. F. Muzy, Eur. 



Phys. J. B 1, 129 (1998); (b) Y. Malecot et al, Eur. 
Phys. J. B 16, 549 (2000); (c) J. Delour, J. F. Muzy, and 
A. Arneodo, Eur. Phys. J. B 23, 243 (2001). 



